package fem3d;

import inf.jlinalg.SolveFailedException;

import java.io.IOException;
import java.util.ArrayList;

import fem2.Constraint;
import fem2.Debugger;
import fem2.Force;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshGenerator;
import fem2.MeshUtilities;
import fem2.Model;
import fem2.Node;
import fem2.analysis.Analysis;
import fem2.analysis.NonlinearStaticAnalysis;
import fem2.element.NonlinearStructuralElement;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.SolverType;
import fem2.enu.State;
import fem2.material.StVenantKirchhoffMM;
import fem2.observer.StrainEnergyObserver;
import fem2.observer.TimingObserver;
import fem2.pre_and_post.GidMeshGenerator;
import fem2.pre_and_post.GidPostStructural;

/**
 * desmonstrate 10-node tetrahedra capability with geometrical nonlinearity and
 * St-Venant Kirchhoff material
 * 
 * @author hbui
 * 
 */
public class TriPodNonlinear {
	static String projDir = "/home/hbui/kratos_janosch";
	static String projName = "tripodnonlinear.gid";
	String fn = projDir + '/' + projName + '/' + "tripod.GiD.msh";

	public Model createModel() {
		MeshGenerator mg = new GidMeshGenerator(fn);
		Mesh mesh = mg.getMesh(3);

		Model m = new Model(mesh);

		mesh.setNumberOfFields(3);

		/*
		 * element creation
		 */
		double E = 1e3;
		double nu = 0.3;
		double t = 1.0;
		State ss = State.THREE_DIMENSIONAL;
		MaterialModel mm = new StVenantKirchhoffMM(E, nu, t, 0, ss);
		for (int i = 0; i < mesh.countBlocks(); i++) {
			StructuralElement e = new NonlinearStructuralElement(mesh.getBlock(i), mm);
			m.addElement(e);
		}

		// mesh.Print();

		/*
		 * conditions
		 */
		Constraint c1 = new Constraint(false, false, false);
		// Constraint c2 = new Constraint(true, false, true);
		ArrayList<Node> bottom = MeshUtilities.seekNodesOnSurface(mesh, 0.0, 1.0, 0.0, 10.0);
		mesh.setConstraint(c1, bottom);
		// mesh.getNode(0).setConstraint(c1);
		// mesh.getNode(1).setConstraint(c1);
		MeshUtilities.correctMesh(mesh);

		double F = -10.0;
		ArrayList<Node> top = MeshUtilities.seekNodesOnSurface(mesh, 0.0, 1.0, 0.0, -13.0);
		for (int i = 0; i < top.size(); i++) {
			Force f = new Force(top.get(i), 0.0, F, 0.0);
			m.addLoad(f);
		}

		return m;
	}

	public static void main(String[] args) throws SolveFailedException, IOException {
		Model M = new TriPodNonlinear().createModel();

		// Mesh mesh = M.getMesh();
		// // Solver solver = new BaitschSolver();
		// Solver solver = new UmfPackLSESolver();
		// solver.setEchoLevel(EchoLevel.OUTPUT);
		//
		// M.preAnalysis();
		// M.Initialize();
		// int n = mesh.countDofs();
		//
		// solver.setSize(n);
		// SymbolicTable t = new SymbolicTable(n, n);
		// t.initializeFrom(M);
		// solver.setSize(t);
		//
		// M.assembleKt(solver.getA());
		// double[] r = new double[n];
		// r = M.assembleR(r);
		//
		// System.out.println("norm r = " + BLAM.norm(BLAM.NORM_TWO, new
		// ArrayVector(r)));
		//
		// double[] u = r.clone();
		// solver.solve(u);
		// M.setX(1, u);
		// M.FinalizeNonlinearIteration();

		/*
		 * analysis
		 */
		// M.setAssembler(new ParallelAssembler(M, 4));
		M.setEchoLevel(EchoLevelType.AUX1);
		Analysis an = new NonlinearStaticAnalysis(M, 2);
		// ((NewtonRaphsonBasedStrategy) an.getStrategy())
		// .setMaxNewtonRaphsonIteration(1);
		// an.getStrategy().setSolverType(SolverType.UMFPACK_SOLVER);
		an.getStrategy().setSolverType(SolverType.MKL_SOLVER);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.AUX1);
		an.getStrategy().addObserver(new TimingObserver());
		an.getStrategy().addObserver(new StrainEnergyObserver(M));
		an.run();

		/*
		 * post process
		 */
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// System.out.println(ArrayFormat.format(mesh.getNode(i).getUHatLocal()));
		// }

		GidPostStructural post = new GidPostStructural(projDir, projName);
		post.clearAllPostFiles();
		post.writeMesh(M, 0);
		post.writeNodalDisplacementResults(M, 0);

		Debugger.warn("Analysis completed");
	}
}
